library(dplyr)
library(tidyverse)
library(padr)
library(readstata13)
library(scales)
library(survey) 

# Specify the directory where the plots should be saved

setwd("insert-your-directory-here")

##############################
# Replicating top-left panel #
##############################

# USA

cc_gal <- read.dta13("Input_Graphs_Gallup.dta", missing.type = TRUE) 

# Europe

cc_eb <- read.dta13("Input_Graphs_Eurobarometer.dta") 

# Create the graph 

pdf("Figure_A1_Panel_A.pdf") 

plot(cc_gal$year, cc_gal$mean_serious_wght, type = "n", xlim=c(2001, 2019),
     ylim=c(0.4, 1), yaxt="n", xaxp=c(2001, 2019,6), 
     xlab="Year", ylab="Share Concerned About Climate Change"
) 

yticks_val <- c(0.4, 0.6, 0.8, 1)

axis(2, at=yticks_val, lab=percent(yticks_val))

# Lines USA

lines(cc_gal$year, cc_gal$mean_serious_wght, lwd=3, col="grey20")

# Lines Europe

lines(cc_eb$year, cc_eb$mean_serious_wght, lwd=3, col="grey20", lty=2)

# Arrows

text(2011.25, 0.81, "Europe", col=1, cex=1)

text(2011, 0.60, "US", col=1, cex=1)

dev.off()

###############################
# Replicating top-right panel #
###############################

europe <- read.dta13("Input_Elections_Europe.dta") 

cmp <- read.dta13("Input_CMP.dta") 

cmp<-cmp[!(cmp$year<"1998"),]

cmp <- subset(cmp, country==42 | country==21 | country==13 | country==14 | country==31
              | country==41 | country==34 | country==53 | country==51 | country==35
              | country==33 | country==32 | country==22 | country==12 | country==11)

# Create the function for the weighted mean

weighted.median <- function(x, w, na.rm=TRUE, ties="mean") {
  if (missing(w))
    w <- rep(1, length(x));
  
  # Remove values that are NA's
  if (na.rm == TRUE) {
    keep <- !(is.na(x) | is.na(w));
    x <- x[keep];
    w <- w[keep];
  } else if (any(is.na(x)))
    return(NA);
  
  # Assert that the weights are all non-negative.
  if (any(w < 0))
   stop("Some of the weights are negative; one can only have positive
        weights.");
  
  # Remove values with weight zero. This will:
  #  1) take care of the case when all weights are zero,
  #  2) make sure that possible tied values are next to each others, and
  #  3) it will most likely speed up the sorting.
  n <- length(w);
  keep <- (w > 0);
  nkeep <- sum(keep);
  if (nkeep < n) {
    x <- x[keep];
    w <- w[keep];
    n <- nkeep;
  }
  
  # Order the values and order the weights accordingly
  ord <- order(x);
  x <- x[ord];
  w <- w[ord];
  
  wcum <- cumsum(w);
  wsum <- wcum[n];
  wmid <- wsum / 2;
  
  # Find the position where the sum of the weights of the elements such that
  # x[i] < x[k] is less or equal than half the sum of all weights.
  # (these two lines could probably be optimized for speed).
  lows <- (wcum <= wmid);
  k  <- sum(lows);
  
  # Two special cases where all the weight are at the first or the
  # last value:
  if (k == 0) return(x[1]);
  if (k == n) return(x[n]);
  
  # At this point we know that:
  #  1) at most half the total weight is in the set x[1:k],
  #  2) that the set x[(k+2):n] contains less than half the total weight
  # The question is whether x[(k+1):n] contains *more* than
  # half the total weight (try x=c(1,2,3), w=c(1,1,1)). If it is then
  # we can be sure that x[k+1] is the weighted median we are looking
  # for, otherwise it is any function of x[k:(k+1)].
  
  wlow  <- wcum[k];    # the weight of x[1:k]
  whigh <- wsum - wlow;  # the weight of x[(k+1):n]
  if (whigh > wmid)
    return(x[k+1]);
  
  if (is.null(ties) || ties == "weighted") {  # Default!
    (wlow*x[k] + whigh*x[k+1]) / wsum;
  } else if (ties == "max") {
    x[k+1];
  } else if (ties == "min") {
    x[k];
  } else if (ties == "mean") {
    (x[k]+x[k+1])/2;
  } else if (ties == "both") {
    c(x[k], x[k+1]);
  }
}


temp <- cmp %>% drop_na(pervote)

temp$year2 <- round(temp$year/5,0)*5

cog_cntry <- tapply(seq(along=temp$environmentalism), list(temp$country, temp$year2),
                    
                    function(i, x=temp$environmentalism, w=temp$pervote/100)  
                      
                      weighted.mean(x[i], w[i]))

cog = colMeans(cog_cntry[,],na.rm = TRUE)

cog <- cog[as.numeric(names(cog))<2021]

cog <-cog / cog[[1]][1]

names(cog)[names(cog) == "2020"] <- "2019"

europe<-europe[!(europe$year<"1998"),]

country_green <- europe %>% group_by(year, country) %>% summarise(country_green = mean(green_share))

country_green$year2 <- round(country_green$year/5,0)*5

green <- country_green %>% group_by(year2) %>% summarise(green = mean(country_green))

bind <- cbind(green, cog)

myvars <- c("year2", "green", "cog")
bind <- bind[myvars]

bind$year2[bind$year2 == "2020"] <- "2019"

Senate <- read.dta13("Input_Senate.dta") 

Senate<-Senate[!(Senate$year<"1998"),]

Senate<-Senate[!(Senate$year=="2017"),]

Senate$year2 <- round(Senate$year/5,0)*5

Senate <- Senate %>% group_by(year2) %>% summarize(senate = mean(COG, na.rm = TRUE))

Senate$senate <-Senate$senate /  Senate$senate [[1]][1]

Senate$year2[Senate$year2 == "2020"] <- "2018"

# Create the graph 

pdf("Figure_A1_Panel_B.pdf") 

par(mar=c(5, 4, 4, 6) + 0.1)

plot(bind$cog, bind$year2,type="n", xlim=c(2000, 2019),
     
     ylim=c(0.5, 1.5),
     xlab="Year", ylab="Environmentalism Indices"
)

axis(side = 1, at = c(2019))

# Add Environmentalism : black dashed line

lines(bind$year2, bind$cog, lwd=4,  lty=2, col="black")  

text(2012.7, 0.84, "Europe: Environmentalism", col=1, cex=1)

# Add Senate : black solid line

lines(Senate$year2, Senate$senate,lwd=3, col="black")

text(2012.5, 1.22, "US: Environmentalism", col=1, cex=1)

# Add second y axis 

par(new = TRUE)
plot(bind$year2,bind$green, type = "n", axes = FALSE,
     ylim=c(0,0.1),  
     yaxt="n",  xlab = "", ylab = "" ) 

yticks_val <- c(0.00, 0.02, 0.04, 0.06, 0.08, 0.1) 

axis(side=4,
     col.axis="red", at=yticks_val,lab =percent(yticks_val))
mtext("Green Share", side=4, col="red", line=3)

# Add Green Share: red dashed line

lines(bind$year2,bind$green,lwd=4, lty=2, col="red")  

text(2004.1, 0.045, "Europe: Green Share", col="red", cex=1)

dev.off()

############################
# Replicating bottom panel #
############################

# Import the dataset

trade <- read.dta13("Input_GDP_Imports.dta")

# Obtain imports over GDP in USA and EU

trade$imports_gdp_us <-trade$imports_us / trade$gdp_usa

trade$imports_gdp_eu <-trade$imports_europe / trade$gdp_eu 

# Set 1995 = 1

trade$imports_gdp_us <-trade$imports_gdp_us / trade$imports_gdp_us [[1]][1] 

trade$imports_gdp_eu <-trade$imports_gdp_eu / trade$imports_gdp_eu [[1]][1]  

# Plot 

pdf("Figure_A1_Panel_C.pdf") 

plot(trade$year, trade$imports_gdp_us, type = "n", xlim=c(1995, 2019),
     ylim=c(1,1.8), yaxt="n", xaxp=c(1995, 2019, 4), 
     main="",
     xlab="Year", ylab="Imports over GDP"
) 

yticks_val <- c(1,1.2,1.4,1.6,1.8)
xticks_val <- c(1995,2001,2007,2013,2019)

axis(2, at=yticks_val, lab=(yticks_val))

lines(trade$year, trade$imports_gdp_eu, lwd=3, col="grey20", lty=2)

lines(trade$year, trade$imports_gdp_us, lwd=3, col="grey20", lty=1)

text(2014.7, 1.5, "Europe", col="black", cex=1)

text(2016.61, 1.3, "US", col="black", cex=1)

dev.off()
